

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h>

#define alpha 1
#define beta 1
double A[N][N +20];
double B[N][N +20];
double x[N];
double u1[N];
double u2[N];
double v1[N];
double v2[N];
double w[N];
double y[N];
double z[N];

void init_arrays()
{
    int i, j;
    for (i=0; i<N; i++) {
        u1[i] = i;
        u2[i] = (i+1)/N/2.0;
        v1[i] = (i+1)/N/4.0;
        v2[i] = (i+1)/N/6.0;
        y[i] = (i+1)/N/8.0;
        z[i] = (i+1)/N/9.0;
        x[i] = 0.0;
        w[i] = 0.0;
        for (j=0; j<N; j++) {
            A[i][j] = ((double) i*j)/N;
        }
    }
}

double rtclock()
{
  struct timezone tzp;
  struct timeval tp;
  int stat;
  gettimeofday (&tp, &tzp);
  return (tp.tv_sec + tp.tv_usec*1.0e-6);
}

int main()
{
  init_arrays();

  double annot_t_start=0, annot_t_end=0, annot_t_total=0;
  int annot_i;

  for (annot_i=0; annot_i<REPS; annot_i++)
  {
    annot_t_start = rtclock();



                                    

/*@ begin PerfTuning (
  def build
  {
  arg command = 'icc';
  arg options = '-fast -openmp -I/usr/local/icc/include -lm';
  }

  def performance_counter
  {
  arg method = 'basic timer';
  arg repetitions = 4;
  }  
  
  def performance_params
  {  
    param T1_1[] = [512];
    param T1_2[] = [1];
    param T2_1[] = [8];
    param T2_2[] = [4];

    param U1[] = [1];
    param U2[] = [4];
    param U3[] = [4];
    param U4[] = [1];

    param PERM1[] = [
      [0,1],
#      [1,0],
    ];
    param PERM2[] = [
     ['i','j'],
#     ['j','i'],
    ];

    param PAR[] = [True];
    param SCREP1[] = [True];
    param SCREP2[] = [True];
    param IVEC1[] = [True];
    param IVEC2[] = [True];
  }
  
  def search
  {
  arg algorithm = 'Exhaustive';
# arg algorithm = 'Simplex';
# arg algorithm = 'Random';
# arg time_limit = 10;
  arg total_runs = 1;
  }

  def input_params
  {
  param N[] = [10000];
  param alpha = 1;
  param beta = 1;

  decl in static double A[N][N+20] = 0;
  decl in static double B[N][N+20] = 0;
  decl in static double x[N] = 0;
  decl in static double u1[N] = 0;
  decl in static double u2[N] = 0;
  decl in static double v1[N] = 0;
  decl in static double v2[N] = 0;
  decl out static double w[N] = 0;
  decl in static double y[N] = 0;
  decl in static double z[N] = 0;
  }
) @*/

/**-- (Generated by ancc)
Input problem sizes:
  N = 10000 
  alpha = 1 
  beta = 1 
Best performance parameters:
  IVEC1 = True 
  IVEC2 = True 
  PAR = True 
  PERM1 = [0, 1] 
  PERM2 = ['i', 'j'] 
  SCREP1 = True 
  SCREP2 = True 
  T1_1 = 512 
  T1_2 = 1 
  T2_1 = 8 
  T2_2 = 4 
  U1 = 1 
  U2 = 4 
  U3 = 4 
  U4 = 1 
--**/

   

register int i,j,k,t;  
register int it,jt,kt,tt;  
register int c1t, c2t, c3t, c4t, c5t, c6t, c7t, c8t, c9t, c10t, c11t, c12t;  
register int newlb_c1, newlb_c2, newlb_c3, newlb_c4, newlb_c5, newlb_c6,  
  newlb_c7, newlb_c8, newlb_c9, newlb_c10, newlb_c11, newlb_c12;  
register int newub_c1, newub_c2, newub_c3, newub_c4, newub_c5, newub_c6,  
  newub_c7, newub_c8, newub_c9, newub_c10, newub_c11, newub_c12;  

/*@ begin PolySyn(  
  parallel = PAR; 
  tiles = [T1_1,T1_2,T2_1,T2_2];
  permut = PERM1;
  unroll_factors = [U1,U2]; 
  scalar_replace = SCREP1; 
  vectorize = IVEC1; 
 
  profiling_code = 'gemver_profiling.c'; 
  compile_cmd = 'gcc'; 
  compile_opts = '-lm'; 
  ) @*/ 



#define ceild(n,d)  ceil(((double)(n))/((double)(d)))
#define floord(n,d) floor(((double)(n))/((double)(d)))
#define max(x,y)    ((x) > (y)? (x) : (y))
#define min(x,y)    ((x) < (y)? (x) : (y))

		
	int c1, c2, c3, c4, c5, c6, c7;

	register int lb, ub, lb1, ub1, lb2, ub2;
/* Generated from PLuTo-produced CLooG file by CLooG v0.14.1 64 bits in 0.02s. */
	lb1=0;
	ub1=floord(N-1,4096);
#pragma omp parallel for shared(lb1,ub1) private(c1,c2,c3,c4,c5,c6,c7)
	for (c1=lb1; c1<=ub1; c1++) {
/*@ begin Loop(
transform Composite(
permut = [['c4', 'c5']],
  regtile = (['c4', 'c5'],[1, 4]),
  scalarreplace = (True, 'double'),
  vector = (True, ['ivdep','vector always']))

  for (c2=0;c2<=floord(N-1,4);c2++) {
    for (c3=max(8*c1,0);c3<=min(8*c1+7,floord(N-1,512));c3++) {
      for (c4=max(4*c2,0);c4<=min(N-1,4*c2+3);c4++) {
        for (c5=max(512*c3,0);c5<=min(N-1,512*c3+511);c5++) {
          B[c4][c5]=u2[c4]*v2[c5]+u1[c4]*v1[c5]+A[c4][c5] ;
          x[c5]=beta*B[c4][c5]*y[c4]+x[c5] ;
        }
      }
    }
  }

) @*/{
  for (c2=0; c2<=floord(N-1,4); c2++ ) {
    for (c3=max(8*c1,0); c3<=min(8*c1+7,floord(N-1,512)); c3++ ) {
      for (c4=max(4*c2,0); c4<=min(N-1,4*c2+3); c4++ ) {
        double scv_11, scv_12, scv_13;
        scv_11=u2[c4];
        scv_12=u1[c4];
        scv_13=y[c4];
        register int cbv_1, cbv_2;
        cbv_1=max(512*c3,0);
        cbv_2=min(N-1,512*c3+511)-3;
#pragma ivdep
#pragma vector always
        for (c5t=cbv_1; c5t<=cbv_2; c5t=c5t+4) {
          double scv_1, scv_2, scv_3, scv_4, scv_5, scv_6, scv_7, scv_8;
          scv_1=x[(c5t+3)];
          scv_2=B[c4][c5t];
          scv_3=B[c4][(c5t+3)];
          scv_4=B[c4][(c5t+2)];
          scv_5=x[c5t];
          scv_6=B[c4][(c5t+1)];
          scv_7=x[(c5t+1)];
          scv_8=x[(c5t+2)];
          scv_2=scv_11*v2[c5t]+scv_12*v1[c5t]+A[c4][c5t];
          scv_6=scv_11*v2[(c5t+1)]+scv_12*v1[(c5t+1)]+A[c4][(c5t+1)];
          scv_4=scv_11*v2[(c5t+2)]+scv_12*v1[(c5t+2)]+A[c4][(c5t+2)];
          scv_3=scv_11*v2[(c5t+3)]+scv_12*v1[(c5t+3)]+A[c4][(c5t+3)];
          scv_5=beta*scv_2*scv_13+scv_5;
          scv_7=beta*scv_6*scv_13+scv_7;
          scv_8=beta*scv_4*scv_13+scv_8;
          scv_1=beta*scv_3*scv_13+scv_1;
          x[(c5t+3)]=scv_1;
          B[c4][c5t]=scv_2;
          B[c4][(c5t+3)]=scv_3;
          B[c4][(c5t+2)]=scv_4;
          x[c5t]=scv_5;
          B[c4][(c5t+1)]=scv_6;
          x[(c5t+1)]=scv_7;
          x[(c5t+2)]=scv_8;
        }
        register int cbv_3;
        cbv_3=min(N-1,512*c3+511);
#pragma ivdep
#pragma vector always
        for (c5=c5t; c5<=cbv_3; c5=c5+1) {
          double scv_9, scv_10;
          scv_9=x[c5];
          scv_10=B[c4][c5];
          scv_10=scv_11*v2[c5]+scv_12*v1[c5]+A[c4][c5];
          scv_9=beta*scv_10*scv_13+scv_9;
          x[c5]=scv_9;
          B[c4][c5]=scv_10;
        }
      }
    }
  }
}
/*@ end @*/

}
/* End of CLooG code */

/*@ end @*/

for (i=0; i<=N-1; i++)
  x[i] = x[i] + z[i];

/*@ begin Loop(
 transform Composite(
  permut = [PERM2],
  regtile = [['i','j'],[U3,U4]],
  scalarreplace = (SCREP2, 'double', 'itv_'),
  vector = (IVEC2, ['ivdep','vector always'])
 )
 for (i=0; i<=N-1; i++)
   for (j=0; j<=N-1; j++)
     w[i] = w[i] + alpha* B[i][j]*x[j];
     ) @*/{
  for (it=0; it<=N-4; it=it+4) {
    double itv_2, itv_3, itv_4, itv_5;
    itv_2=w[(it+2)];
    itv_3=w[(it+1)];
    itv_4=w[(it+3)];
    itv_5=w[it];
    register int cbv_1;
    cbv_1=N-1;
#pragma ivdep
#pragma vector always
    for (j=0; j<=cbv_1; j++ ) {
      double itv_1;
      itv_1=x[j];
      itv_5=itv_5+alpha*B[it][j]*itv_1;
      itv_3=itv_3+alpha*B[(it+1)][j]*itv_1;
      itv_2=itv_2+alpha*B[(it+2)][j]*itv_1;
      itv_4=itv_4+alpha*B[(it+3)][j]*itv_1;
    }
    w[(it+2)]=itv_2;
    w[(it+1)]=itv_3;
    w[(it+3)]=itv_4;
    w[it]=itv_5;
  }
  for (i=it; i<=N-1; i=i+1) {
    double itv_6;
    itv_6=w[i];
    register int cbv_2;
    cbv_2=N-1;
#pragma ivdep
#pragma vector always
    for (j=0; j<=cbv_2; j++ ) {
      itv_6=itv_6+alpha*B[i][j]*x[j];
    }
    w[i]=itv_6;
  }
}
/*@ end @*/

/*@ end @*/

    annot_t_end = rtclock();
    annot_t_total += annot_t_end - annot_t_start;
  }
  
  annot_t_total = annot_t_total / REPS;
  printf("%f\n", annot_t_total);

  return ((int) w[0]); 

}
